Dynamics of parametric matter wave amplification 
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We develop a model for parametric amplification, based on a density matrix approach, which 
naturally accounts for the peculiarities arising for matter waves: significant depletion and explicit 
time-dependence of the source state population, long interaction times, and spatial dynamics of the 
amplified modes. We apply our model to explain the details in an experimental study on twin-atom 
beam emission from a one-dimensional degenerate Bose gas. 
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I. INTRODUCTION 

Parametric amplification of quantum mode popula- 
tions is a feature common to both light, and matter 
waves. For light, optical parametric amplification [TJ[2] is 
the key technique to populate twin-modes containing cor- 
related photon pairs. For atomic matter waves, numerous 
emission and amplification schemes have been demon- 
strated both in spontaneous f^-'F and Bose-enhanced pl- 
[T2] regimes. Ideally, the populations of two modes, se- 
lected by fundamental conservation laws, grow identically 
during the amplification process, analogous to the signal 
and idler modes in parametric down-conversion of pho- 
tons. Key features of such twin-atom beams, such as sup- 
pressed relative and enhanced absolute number fluctua- 
tions l9Ul2| , non-trivial second-order correlations [H |8] 
and quadrature squeezing |10fll2| have been shown ex- 
perimentally. 

A crucial difference between matter and photon twin- 
beams arises from the microscopic process driving the 
stimulated emission. As photons do not interact, the am- 
plification process has to be mediated by a medium which 
is being pumped by a strong light field. Due to the rela- 
tively weak X nonlinearity and the short, localized in- 
teraction in the medium, the conversion efficiency is low, 
allowing to neglect the depletion of the pump beam and 
to reduce the description to the signal and idler modes 
only (undepleted pump approximation). For interacting 
matter waves, the pump field itself acts as the nonlin- 
ear medium. Interaction times can be long, so that the 
depletion of the source state considerably affects the dy- 
namics. Furthermore, the coupling between source and 
amplified modes typically extends over the entire sys- 
tem size, and the mode structure can be more complex. 
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At finite temperatures, thermal phase fluctuations reduce 
the coherence, effectively depleting the source state [13]. 
Finally, the source population may explicitly depend on 
time, e.g. if the emission process starts before pumping is 
completed, the source is being replenished continuously, 
or other loss channels are present. 

In this article, we develop a simple and numerically 
tractable model for the description of stimulated matter 
wave emission into twin-beams. Our approach is based 
on density matrices |14| , and extends previous theoretical 
studies which have concentrated on the regime of rapid 
scattering into many, weakly occupied modes, where the 
Bogoliubov approximation (equivalent to negligible de- 
pletion) holds [15 -17]. In contrast to calculations based 
on the positive-P [ISfETj or truncated Wigner [22j meth- 
ods, our approach does not rely on stochastic sampling, 
and remains valid for long interaction times and arbitrary 
mode populations. We will be primarily concerned with 
modeling of de-excitation experiments from a vibrational 
state [31, where we record the population of amplified 
twin-beams over time, but our density matrix calcula- 
tion can be applied to a much larger class of twin-beam 
experiments. 
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FIG. 1. (Color online) Overview of the twin-atom process. 
(a) Level scheme of the involved states and transitions, (b) 
Qualitative population dynamics of the states as shown in (a) 
once the pumping process starts. 
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FIG. 2. (Color online) Population transients Pg,p + Pg,-p, Pc.o 
for the two-mode model of eq. M. In the simulations we use 
700 atoms. Thick and thin lines report the normalized popu- 
lation of \g, ±p) and |e, 0), respectively. Solid and dashed lines 
correspond to simulations performed for the full Schrodinger 
equation and the density matrix approximation of eqs. (|3]|4|. 
The dotted line indicates the growth of the emitted popula- 
tion for an undepleted pump. 



II. THEORY 

We consider a trapped, degenerate Bose gas in one 
dimension, such as the ground state \g, 0) of a very elon- 
gated harmonic trap with weak confinement Wj. along 
its longitudinal direction x and negligible occupation of 
transversally excited states. Starting from t = 0, we 
pump the system into a source state |e, 0), carrying a 
per-particle excess energy e ^ huj^ but leaving the spa- 
tial wave function along x unchanged (fig. IT]). This energy 
can be stored in internal degrees of freedom [15], dimer 
molecules [18] . or in a vibrational state orthogonal to 
X [HI [53] as in the experiment described below. Conserva- 
tion laws for parity and momentum dictate that the state 
|e, 0) can only decay via a two-body process into paired, 
propagating modes \g, ±p) with identical population and 
opposite momenta centered around zLp « ±v 2?ne, where 
m denotes the particle mass. A slight shift may arise 
from mean-field effects |53]. For a scheme populating 
twin- modes in a spinor gas [lOflT^ . the equivalent states 
would be the |m = 0) state as a source and two states 
with opposite magnetizations |±r7i) instead of \g,±p). 



A. Two-mode model 

If we consider a single pair of modes resonant with |e, 0) 
and instantaneous excitation at i = 0, the relevant part 
of the two-body interaction Hamiltonian can be written 
as 
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(1) 

similar to parametric amplification in optics, where the 
strength of atom-atom interaction k corresponds to the 



nonlinear susceptibility. 

A full numerical solution of the Schrodinger equation 
with the two-mode Hamiltonian (IT]) can be accomplished 
by expanding the many-body wave function in terms of 
bosonic Fock states, where m atom pairs are promoted 
from |e, 0) to \g,±p): 
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Here Cm are the wave function amplitudes in Fock space 
and |0) is the vacuum state. The Schrodinger equation 
with the Hamiltonian H^m is solved with the Crank- 
Nicolson method according to the prescription given in 
ref. [25_. 

Solid lines in fig. [2] show results for a simulation of 
n = 700 atoms, which initially all reside in the transver- 
sally excited |e, 0) state. The non-linear interaction of 
eq. (IT]) promotes atoms pairwise to the twin-atom states 
|(7, ±p). Initially the process is slow and governed by 
spontaneous scatterings. Only when a sufficient popula- 
tion has built up in \g, ±p), say at times around 2/ nn, the 
scattered atoms act as a seed for the ensuing rapid stim- 
ulated emission which continues until the |e,0) reservoir 
is emptied. Finally, owing to the second term in paren- 
theses of eq. ([T]) , atoms scatter back from the twin-atom 
states to |e,0). 

Unfortunately, for many near-resonant modes, the ex- 
ponentially increasing size of the Hilbert space prohibits 
a generalization of this direct approach. To overcome 
this problem, we introduce a density-matrix description 
for the two-level problem, which we will extend to many 
modes in the next section. The lowest moments are the 
densities pe,o = (al,o"e,o) ^^^ pg^p = (ajpflg^). Their 
dynamic equations can be obtained from the Heisenberg 
equations of motion, 

Pg,±p = -^Pefl = 2k 3(A) , A = (4.p4.-p(^e.o) ) ' 

(3) 
Here, the densities are driven by the two-particle coher- 
ence A, whose time evolution is in turn governed by 
three-particle couplings. To truncate this hierarchy of 
equations of motion, we introduce a correlation expan- 
sion in the spirit of ref. [2], and factorize all three- 
particle density matrices into densities of lower order. 
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Initially, pg^p — and peo — n ^ 1, and eq. (W| re- 
duces to A sa iKn^. Consequently, the emitted popu- 
lation grows quadratically as Pg,±p ~ {nnt)'^, which for 
short times t is consistent with the result for an unde- 
pleted pump mi nil: p^g^p = sinh^(Kni). In eq. (]4|, 
the second term in brackets accounts for the backscat- 
tering of population into the source state. A comparison 



of the results of the density matrix equations (131 m , see 
dashed lines in fig. |2] with the exact solution shows good 
agreement. 



B. Multi-mode model 

In a realistic matter-wave amplifier, several twin- 
modes become populated simultaneously, as the amph- 
fication bandwidth is broadened by the mean field of the 
source. In contrast to a full wave function approach, it is 
numerically straightforward to extend the density matrix 
framework of eqs. (3]4| to a multi-mode Hamiltonian 
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where k^j is the interaction matrix element between the 
source and the different spatial modes \g,i). Momentum 
and parity conservation are fulfilled due to the symme- 



tries of the matrix elements 



The number of neces- 



sary modes M can be estimated from the per-particle 
mean field energy /z of the source [16, : AI '^ fi/hto, where 
huj is the typical energy level spacing of adjacent modes. 
The multi-mode description inherently includes spatial 
dynamics, such as the propagation of twin-beam wave 
packets. 

As for the two-mode case outlined in the previous sec- 
tion, we employ a factorization scheme for many-particle 
density matrices to obtain the density-matrix equations 
for the many-mode Hamiltonian ([5]). Let Pg,y be the 
single-particle reduced density matrix for the twin-atom 
states, generalizing Pg^±p of the two-mode case. The dy- 
namic equations for p directly follow from the Heisenberg 
equations of motion 
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Here A^j = (a^aUao)^) is the two-particle coherence be- 
tween source and emitted modes, whose time evolution 
can be evaluated to 
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With this result, the dynamics of the twin-mode popula- 
tions can be derived, similarly to the two-mode case. 

III. EXPERIMENT 

As a representative example of twin-beam creation, in 
this paper we consider a coUisional de-excitation scheme 




FIG. 3. (Color online) Population of twin-atom modes, nor- 
malized to the total atom number. The data points with 
error bars report the experimental results for optimal (IV) 
and scaled excitation (other lines) with scaling factors of ap- 
proximately 0.3 (I), 0.5 (II), 0.8 (III), and 1.2 (V). Solid lines: 
corresponding theory results X^fcPs.fefe- Dashed line: theory 
for optimal excitation, where additional spontaneous emis- 
sion processes are neglected. Dash-dotted line: Upper bound 
imposed by the population of the Onsager- Penrose mode. (In- 
set) Variance of emitted fraction, relative to shot noise (sn) , 
for optimal excitation (solid, black) and scaling 0.3 (dotted, 
blue). Each point is an average over seven adjacent times 
and corrected for imaging noise |26| and total atom number 
fluctuations. 



based on a one-dimensional, quantum degenerate Bose 
gas of Rubidium-87 atoms trapped on a chip |27l [25] . 
as described in more detail in ref. [S]. Within 5 ms, 
atoms are pumped into the source state |e,0), which is 
a vibrationally excited transversal state of the confining 
waveguide potential. This pumping is accomplished via 
a mechanical optimal control protocol, populating |e, 0) 
with almost unity efficiency by non-adiabatic translation 
of a strongly anharmonic trapping potential along the 
transversal direction y (see appendix [A|. In fig.^we show 
the population of the twin-beam modes, relative to the 
total atom number n « 800, as a function of time after 
starting the excitation sequence. Experimental points are 
obtained by counting photons in appropriately defined 
regions within fluorescence images taken after 46 ms of 
free expansion, see appendix [B] for details. One observes 
that for the optimal ramp (black markers, series IV) the 
twin-atom population increases over approximately 10 ms 
and finally reaches a plateau value. About 40 percent of 
the atoms are emitted into twin beams. We repeated 
the measurement several times, scaling the amplitude of 
the optimal trap motion used for pumping by different 
factors [29_. For these non-optimal excitation protocols 
(series I-III and V) the final twin-atom populations are 
reduced. 



To model these experiments, we need to describe the 
time-dependent excitation (pumping) of atoms into the 
source state with energy e, see fig. Iljb), taking place on 
a time scale on the order of ^ lOhJe « 5 ms with h the 
Planck constant, which is similar to that of the emission. 
Additionally, we have to account for the thermal deple- 
tion of the source mode and the spatial dynamics of the 
twin beams. As will be demonstrated next, all of this can 
be conveniently accomplished within the density matrix 
approach. 



A. Pumping 

Fig. Qa) shows the time-dependent position of the trap 
minimum (solid line) together with the time-dependent 
density distribution along y, which is obtained from the 
solution of the one-dimensional Gross-Pit aevskii equa- 
tion (GPE) along y. We can approximate the excita- 
tion dynamics by a two-mode description comprising the 
states \g,0) and |e,0), which are chosen as displaced 
eigenstates [30] of the harmonic part of the waveguide po- 
tential with a common, periodic displacement (see panel 




4 6 

Time (ms) 

FIG. 4. (Color online) Numerical twin-beam dynamics re- 
sults, from the simulations for our experiment, (a) Density 
distribution along the transversal direction y. The solid line 
indicates the driven motion of the trap minimum. The dashed 
line shows the collective oscillation frame used to define the vi- 
brational states, (b) Wave functions of the vibrational states 
along y, corresponding to \g,0) ,\g,iizp) and |e,0), respec- 
tively, at zero collective displacement, (c) Populations of \g, 0) 
and |e,0), derived by projecting the numerical results shown 
in (a) and neglecting the twin-beam emission. The dashed 
(green) line shows a constant Rabi coupling model, (d) Illus- 
tration of one particular twin-beam mode ^ip with energy e 
within the effective longitudinal potential, which is defined by 
the trap confinement and the mean field of the source cloud, 
(e) Spatial evolution of the twin-beam modes along the longi- 
tudinal axis, which are created in the center of the condensate 
region and then propagate out. The shaded area corresponds 
to the density distribution of the initial state, the dashed lines 
indicate its radius. 



(b) ). This displacement is reflecting the collective oscil- 
lation (fig.llFa), dashed line), which is necessary to drive 
the system into the strongly anharmonic part of the po- 
tential where the state transfer occurs, but is otherwise 
irrelevant for the dynamics between the states. Project- 
ing the GPE results onto the displaced state basis yields 
the population of \g, 0) and |e, 0), as depicted in panel (c). 
A near-resonant Rabi coupling with an initial time delay 
(dashed line in panel (c) ) can be fitted to the population 
dynamics to provide a suitable model for the following 
calculation steps. The optimal control sequence (lead- 
ing to series IV in fig. Isl then corresponds to a 7r-pulse 
which completely inverts the system from \g,0) to |e, 0). 
See appendix |C] for a more comprehensive discussion of 
the transversal dynamics. 



B. Thermal source depletion 

We now turn to the emission dynamics along the lon- 
gitudinal direction x. Especially for an elongated sys- 
tem, as studied here, thermal longitudinal phase fiuc- 
tuations are expected to have a strong influence on the 
experimental results by effectively depleting the source 
mode [13J. We assume that the atoms are initially in 
the Thomas-Fermi ground state of a harmonic trap along 
X, and account for phase fluctuations through a den- 
sity matrix p{x,x') obtained from quasi-condensate the- 
ory |2I], where the Bogoliubov excitations are populated 
thermally. Next, we split p into a condensate (Onsager- 
Penrose mode) and a thermally excited part pth{x,x'), 



p{x, x') = no(l>Q{x)(j)o{x') + pth{x, x') 



(8) 



where ng is the largest eigenvalue of the density matrix 
and 4'o{x) the corresponding eigenfunction. For our cal- 
culation we assume T = 25 nK, which is compatible with 
experimental observations [32' and leads to no w 0.44 n 
with n fa 800 (see dash-dotted line in flg. [3]). As the 
twin-atom production is driven by a flxed phase relation 
between the source and twin-atom states, eq. pi, the 
non-condensed part is expected to have only little influ- 
ence on the emission dynamics, and will be neglected in 
the following. 



C. Twin-beam states 

To apply the multi-mode description of eq. (|5| to our 
flnite-size system, we compute a set of « 30 highly ex- 
cited (real) single-particle states (pi (x) at energies around 
e ~ h'^kQ/2m, for an effective potential including the 
mean fleld of the source state. One typical state is de- 
picted in flg.kld). The coupling matrix elements of Hmm 
can now be expressed as 



i^ij =9 \'ipgiz)\ dz / ^p*{yf^pe{yfdy 



(9) 



X / (j)*{x)<j)*{x)Mxrdx, 



with i/)e,g(j/) and il'giz) denoting the wave functions of 
the transversal states defined above, and g — Airfi^as/m, 
where a^ and m are scattering length and mass of the 
atoms. The full dynamics of the twin-atom production 
is now governed by the Hamiltonian 



H — Ho + Hp 



H 



MM, 



(10) 



where Hq accounts for the free evolution of |e, 0), \g,0), 
and \g-,i), and ^pump describes the pumping process 
through the Rabi-type excitation described above. 



D. Results: optimal excitation 



For the Hamiltonian in eq. ([10|, we numerically solve 
the equations of motion of the density matrices. Fig. We) 
shows the computed real-space density of twin-atoms as 
a function of time. One observes that the twin-atoms are 
emitted in the condensate region initially, and then prop- 
agate as packets with group velocities -izhko/m, reach- 
ing the edges of the Thomas-Fermi distribution a few 
ms after the end of the excitation pulse. All parameters 
of our model are obtained for realistic trap potentials 
(source excess energy e « ft, x 1.8 kHz) and experimental 
parameters taken from ref. [9 . For the optimal excita- 
tion protocol, the slope and final value of the emission 
dynamics (dashed line in fig. Isla) ) is in good agree- 
ment with the experimental observation (black points), 
with the exception of a time delay of w 2 ms. This in- 
dicates that we underestimate the emission rate at early 
times, where even a very small twin-beam population suf- 
fices to trigger Bose-enhanced emission. Such scattering 
may be caused by the non-condensate source modes, or 
a non-thermal population of high-momentum modes due 
to technical noise causing premature excitation. To ac- 
count for this, we add a weak channel for scattering atoms 
into longitudinal modes with momenta around ±fco which 
act as an additional seed in excess of vacuum fiuctua- 
tions. Good agreement can be reached for a scattering 
rate of F « 0.4s~^ • n^, which is equivalent to typically 
~ 3 • IQ^^Kii ■ rig per mode, with Uf. = {a\ oflco)- Indeed, 
with this addition the simulation (series IV, solid line) 
accurately reproduces the experimental result. 



E. Results: modified excitation 

We next address the results for scaling the amplitude 
of the trap motion by different factors (series I-III, V), 
which again can be well described by Rabi pulses of dif- 
ferent area. The couplings used in the calculations for 



the twin-beam population growth are derived from obser- 
vation of the transversal momentum distribution during 
the excitation (see appendixlcl). For series I, IV, and V in 
fig. [3] this procedure yields good agreement without any 
further free parameters. Only for intermediate excitation 
scalings (series II, HI), a slight adjustment of the exci- 
tation coupling was necessary, which we attribute to the 
approximations made in the pumping model (see above) . 
Note, that for the lowest scaling Bose enhancement is 
weak, which provides a stringent means for determining 
r. The complete set of experimental results can be nicely 
modeled by our theory (solid lines in fig. jsl . 

To illustrate the amplified character of the twin-beam 
creation, we investigate the fluctuations of the twin 
beams population N at different times. In fig. [sl (in- 
set) we show the measured relative variances ^^ ex oj^/ N 
of the relative population of emitted pairs over many ex- 
perimental realizations. They are corrected for imaging 
noise and total atom number fluctuations, and normal- 
ized to an approximation to the shot noise expected for 
random spontaneous emission, neglecting pump deple- 
tion and temperature fluctuations. See appendix [D] for 
details on the variance calculation. For optimal excita- 
tion (black dots, solid line), a pronounced peak near the 
maximum slope of the population growth indicates the 
exponential ampliflcation of initial fluctuations. Such be- 
havior is absent in the experiment with the weakest exci- 
tation (blue dots, dashed line). In contrast, the relative 
number fluctuations of the twin-atom clouds are strongly 
suppressed, as discussed in detail in ref. ^. 



IV. CONCLUSION 

In conclusion, we have derived a density matrix ap- 
proach to quantitatively analyze emission of matter 
waves in the strongly Bose-enhanced regime, which has 
recently been reached in various experiments. Neither 
the Bogoliubov approximation, nor stochastic methods 
are employed, making our approach eligible for strong de- 
pletion of the source and long interaction times. We have 
presented experimental results for twin-beam population 
growth in a one-dimensional degenerate Bose gas, which 
are governed by source depletion, spatial dynamics, and 
explicit time-dependence of the source population. The 
good agreement between experimental results and theory 
predictions suggest that the physics underlying amplifled 
emission of twin-atoms in a real experiment is captured 
by our model. Next steps will comprise further experi- 
mental and theoretical studies on properties of strongly 
populated twin-atom beams beyond single-particle densi- 
ties, with particular regard to second-order correlations, 
revealing number squeezing [5, 9 12 or violation of clas- 
sical inequalities |33| . 
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Appendix A: Trap preparation 

The ultracold Rubiduni-87 gas acting as matter-wave 
source is prepared using laser cooling and then loaded 
into an atom chip wire trap |27l I28| , where forced evap- 
orative cooling to quantum degeneracy occurs. The ini- 
tially transversally symmetric lofFe-Pritchard field con- 
figuration created by the chip wires is modified by radio- 
frequency dressing |34l [55] . Typically used for creating 
double well potentials, this technique also allows for the 
introduction of anharmonicity and anisotropy to a sin- 
gle trap when the dressing strength is kept below the 
point where actual splitting of the potential occurs. Us- 
ing chip wires running parallel to the main trapping wire, 
we apply an ac magnetic field of ^ 0.75 G peak-to-peak 
amplitude at a frequency red-detuned by —54 kHz with 
respect to the atomic Larmor frequency near the trap 
minimum (824 kHz). 

The resulting potential can be calculated numerically 
by means of a Floquet analysis [36 . In the two transver- 
sal directions it can be approximated by quartic poly- 
nomials of the form V/h = ^a;(r/ro)^ -I- A(r/ro)^, where 
To = yjhjrmjj (Duffing oscillator) with the atomic mass 
771. For the y-direction, along which the excitation 
is performed the parameters are ujy /2tt — 1.64 kHz, 
To, J, = 0.266 Jim and \y/2-K — 74.9 Hz. For the z- 
direction perpendicular to the excitation motion the pa- 
rameters are uJz/It: = 2.50 kHz, tq.z = 0.216 jim and 
Xz/'^T^ ~ 25.9 Hz. By solving the Schrodinger equa- 
tion for these potentials we can determine the first trap 
levels along y and z as [Ey^i,Ey^2TEz,i,Ez^2] /"^^t^^ — 
[1.82, 3.78, 2.56, 5.18] kHz, if the respective zero-point en- 
ergies are subtracted. Along the longitudinal x-axis, the 
harmonic trap frequency lo^JIt: = 16.3 Hz is determined 
by observation of a deliberately excited sloshing mode of 
the atom cloud. 

The fast transversal motion of the potential is accom- 
plished by applying a modulated current to an auxiliary 
wire running parallel to the main trapping wire, moving 
the trap minimum by 26 nm/mA along y and 9 nm/mA 
along z. Due to the anisotropy of the transversal po- 
tential, the motion along z is off-resonant and has no 
significant influence. 

The geometry of all involved chip wires and homoge- 
neous offset fields is shown in fig. |5] 



Appendix B: Detection and image analysis 

After starting the excitation sequence, we run the ex- 
citation and emission process for a given time t. At this 
time, we suddenly switch off the trapping potential, im- 
plying that for t < 5 ms the excitation process is still in- 
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FIG. 5. (Color online) Schematic of the atom chip layout. The 
waveguide potential is formed by the current /trap through the 
main trapping wire and a static magnetic field By . On a sep- 
arate chip layer, currents 7ax in broad wires provide axial con- 
finement. An external field B^ completes the loffe-Pritchard 
configuration. The radio frequency dressing currents 7rf are 
applied to wires in parallel to the trapping wire. Finally, the 
modulation of the trap position is accomplished by a current 
/mod in an auxiliary wire. 
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FIG. 6. (Color online) Experimental image analysis. (Inset) 
Typical experimental image data for optimal excitation and 
t < 5.5 ms, averaged over 12 shots. As the image is taken after 
46 ms of expansion time, it predominantly refiects the initial 
momentum distribution. The scale bar corresponds to a spa- 
tial distance of 187 iim, equivalent to the typical momentum 
of 5.5]im~^ ~ ko. The dashed line indicates the transversal 
center-of-mass position of the emitted clouds. (Main plot) 
Longitudinal momentum distribution obtained from integra- 
tion over the image. The shaded area indicates a fitted ther- 
mal distribution with temperature T « 25 nK. Dashed lines: 
range (source region) outside of which atoms are counted as 
emitted population. 



complete. The fast transversal expansion of the cloud due 
to the tight waveguide confinement causes atom interac- 
tions to vanish rapidly, and the ensuing expansion can be 
considered ballistic. After itot — 46 ms of expansion, we 
take a fiuorescence image |26J . fully integrating over the 
z-direction (see inset of fig. l6|. Along y (i.e., integrating 
over x), the resulting image represents the initial mo- 
mentum distribution, as the transversal cloud size before 
expansion is negligible (far field) . If we express momenta 



as wave numbers ky, a distance 6y in the image hence cor- 
responds to Sky = aSy with a = m/httof ~ 0.03 jim"^. 
Along the longitudinal direction x, the far field condi- 
tion is not fully reached due to the initial condensate 
radius of L « 20]jm and the typical momenta in the 
quasi condensate as given by the thermal phase coher- 
ence length l^ fST: k^ = IT^ ~ 0.1 jim"-'^ « aL. Using an 
appropriate model for the initial density and momentum 
distributions (classical fields method |37]), we can fit a 
longitudinal thermal distribution to the expanding cloud 
(fig. pi, yielding an estimate to its temperature. Within 
our current analysis, this relies on the assumption that 
the excitation and emission process does not lead to a 
strong modification of the momentum spectrum at short 
times t. A more detailed study of the interplay between 
thermal effects and the emission dynamics is under way. 
On the other hand, the center momentum of the emit- 
ted atoms is ko sa 5.5 jim"^ 3> {aL,kci,), hence the lon- 
gitudinal overlap of the source cloud and the emitted 
beams is negligible. The fraction of emitted atoms at 
time t (as shown in fig. Isl) can thus be determined by 
simple counting of fluorescence photons in appropriately 
defined ranges of the longitudinal distribution (outside of 
dashed lines in filg.pl). For each combination of t and ex- 
citation amplitude we repeated the experiment typically 
12 times to obtain reliable estimates for the mean and 
variance of the emitted fraction. 



Appendix C: Transversal dynamics 

The vibrationally excited state |e,0) of the anhar- 
monic transversal potential is being populated by os- 
cillatory displacement of its origin V{y — yo{t)), follow- 
ing an optimal control ramp. To determine the func- 
tion yo(t), we propagate the wave function with the one- 
dimensional Gross-Pitaevskii equation (GPE), starting 
from its ground state, and maximize the overlap with 
the first excited state at t = 5 ms. At t > 5 ms, excita- 
tion stops and yo is held constant. As outlined above, 
our experimental time-of-fiight images, integrated over x 
within the bounds shown as dashed lines in fig. [6] di- 
rectly reflect the momentum density distribution along 
y, allowing direct comparison to the results numerically 
obtained from GPE (flg. l7k,b). Excellent agreement can 
be reached, if an independently measured linear response 
function of the electronics employed for the experiments 
is taken into account, which slightly modifles the yo{t) 
trajectory in the experiment. Once the fraction of emit- 
ted twin-beams becomes significant (t ^ 5 ms), the agree- 
ment gradually gets worse, as damping due to the emis- 
sion is not accounted for in the excitation model. For the 
scaled control ramps, we left the amplitude scaling as free 
parameter in the GPE simulation and directly optimized 
agreement with the experimentally observed dynamics. 

As described in the main text, to properly define the 
orbitals of the states \g,0) and |e,0) we transfer into a 
coordinate system that is co-oscillating with the (quasi- 



classical) sloshing of the cloud within the harmonic part 
of the potential V{y). This collective oscillation is in- 
herent to our excitation technique, as it is necessary to 
access the anharmonic part of the potential. Within the 
co-oscillating frame, the anharmonicity acts as an oscil- 
lating force in excess of the restoring force, enabling the 
transfer into an excited eigenstate |38[I39| . Consequently, 
we approximate the wave functions of I5, 0) and |e, 0) 
with the displaced first and second eigenfunctions of a 
harmonic oscillator 1301: 



\g,0)^D[a{t)]\0) 
\e,0)^D[ait)]\l) 
D[a{t)] — exp[a{t)cr — a{t)*a], 



(CI) 



where \n) designates the n-th number state of a harmonic 
oscillator with frequency ujy as defined above, and a, a^ 
the corresponding ladder operators. The complex dis- 
placement a(t) corresponds to the origin y{t),ky(t) of 
the co-oscillating frame in phase space, with 



y ^ V2ro,y^[ait)l 



^^< 



%a(t)\. 



The question remains how to appropriately define the 
co-oscillating frame y{t),ky{t) in our analysis. A strin- 
gent means to obtain its origin is readily provided by ob- 
serving the average transversal momentum of the emitted 
atoms, which are in state \g, ±p) that by construction has 
the same transversal center-of-mass as \g, 0) (see fig.|6] in- 
set), hence fully defining the collective oscillation. How- 
ever, during the first phase of excitation (especially for 
weak excitation scalings) when the emitted atom number 
is low, the uncertainty of the experimentally determined 
center-of-mass is large. Moreover, at later times the 
emission increasingly perturbs the transversal dynamics, 
making a direct extraction of the oscillating frame from 
experimental data inconsistent with our model, which 
omits emission. 

Instead, we proceed by numerically optimizing the 
phase space trajectory a{t) of the co-oscillating frame 
origin, at each time t maximizing the overlap x[a(i)] of 
the GPE result \ilj{t)) with a superposition of the basis 
states: 



X[a{t)]^ m)\D[ait)]\^l-^it)\0) 



with the excited state population ri{t) as free parameter. 
For optimal excitation, this procedure yields x ^ 0-95 
overlap at all times t. In fig. [Tli, a comparison of the nu- 
merically optimized ky and the center-of-mass position 
of the emitted atoms is shown. The excellent agreement 
at times t < 5 ms strongly corroborates the displaced os- 
cillator eigenstate model. In fig. ^, the momentum dis- 
tribution of the two-state superposition is shown, refiect- 
ing the smooth transition from ground to excited state, 
with strong beating between the states at intermediate 
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FIG. 7. (Color online) Transversal momentum dynamics. Left: optimal excitation, right: excitation amplitude scaled by 
« 50 %. (a) Transversal momentum distribution of the excited cloud, obtained from fluorescence images by integrating over 
X within the source region (dashed lines in fig. p|. (b) Numerical result obtained from one-dimensional GPE. For the scaled 
ramp, the scaling amplitude is fitted, whereas there are no free parameters for the optimal ramp. Solid line: numerically 
optimized co-oscillating reference frame position ky{t) as described in the text, (c) Momentum distribution obtained from 
the two-state approximation within the co-oscillating frame. This result corresponds to the effective populations used for the 
emission calculation (fig. ^) . (d) Dots: transversal center-of-mass momentum of the emitted atoms, as extracted from images 
(dashed line in fig. [6] (inset) ). Solid line: same as in (b). 



-I 








IV (opt.) 


0.8 
0.6 


(a) 




if/f 


w 


0.4 




JV 


9, 


- " ^ 


0.2 






y\ 


^' 





^.^ 


ts- 


<r<'-=^ 


\ 



2 3 4 

Time (ms) 




4 6 8 10 
Time (ms) 



FIG. 8. (Color online) (a) Grey lines: populations of the ex- 
cited state r\ extracted from the GPE result within the two- 
state superposition model. Solid lines: fits of a constant Rabi 
coupling model to the numerical result. Dashed lines: popula- 
tions for Rabi coupling parameters that were used as input to 
the emission calculations for intermediate excitation strength 
(green, red), (b) Resulting emission dynamics curves, corre- 
sponding to the excitation curves shown in (a), respectively. 
Experimental data for the intermediate scalings (II and III) 
are shown as dots. 



excited fractions. Fig. Isja) (grey lines) shows the pop- 
ulation of the excited state ?/(i), for the five excitation 
strengths inferred from comparison of the GPE result to 
the experimental data (fig. ItI. The smooth transitions 
observed in the co-oscillating frame are in stark contrast 
to the strong transient population of higher states (as e.g. 
shown in [38 , 39 ) which would be obtained by simple pro- 
jection of the GPE result on the oscillator eigenfunctions. 



To facilitate the inclusion of continuous pumping into 
the emission dynamics calculation, we approximate the 
excitation process by a constant, near-resonant Rabi cou- 
pling with coupling strength f], detuning 8 and initial 
time delay toi ^^nd obtain 



B, 



2 





e,0)(.g,0|+H.c. 



io < t < 5 ms 
otherwise 



pump v^ J 

^ = Q.-cx^{-i{Ey,^lh + S){t-U)\ (C2) 



for the pumping Hamiltonian in eq. (10). This corre 



sponds to an approximation f](t) to the numerically ob- 
tained two-state populations: 



m 



^ [ {§) sin^ [Q.'{t ~ to)] : to < t < 5ms 
() : otherwise 



(C3) 



where VL' — \/fP~+~^. Fig. p[a) (solid lines) shows, that 



eq. ( C3 ) is an appropriate approximation to the numeri- 



cally obtained ri{t), if Jl, A and io are left as fit parame- 
ters. 

For the optimal, weakest and strongest excitations 
(lines IV, I, and V in fig.jsj respectively), using the fitted 
parameters as input to the dynamics calculation directly 
yields excellent agreement of the emitted atom number 
to the experiment as shown in fig. |3] Only for the two 
intermediate strengths (lines II and III), emission curves 
derived in this way show unsatisfactory agreement (solid 
lines II and III in fig. Islb) ). We attribute this to fail- 
ure of the two-level approximation to yield estimates for 
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the ri{t) that are sufficiently precise for prediction of the 
Bose-enhanced emission dynamics, which is highly sen- 
sitive to changes in the initial source population. To be 
able to verify the validity of our emission dynamics the- 
ory for intermediate excitations, only for those two set- 
tings we slightly modified the parameters used as input 
to the emission calculation. The resulting source popu- 
lations fj{t) are shown as dashed lines in fig. p[a). Even 
the slight deviation from the numerically obtained r]{t) is 
sufficient to considerably change the emission dynamics 
results and reach agreement to the experiment (dashed 
lines in fig. [slb) ) . 

A more detailed description of the transversal dy- 
namics would additionally require the explicit treat- 
ment of the orbital dynamics, e.g. within the niulti- 
configurational Hartree method for bosons | .40::41, . which 
is beyond the scope of our present paper. 



Appendix D: Calculation of relative variances 

In fig. [3] (inset), we show the fluctuations of the relative 
population of the twin-beam states in terms of the rela- 
tive variance f^. It can be shown that the fluctuations 
of the emitted number of atoms N can be written with 
respect to the variance of a Poissonian distribution as 



Var(iV - pn) 



2fN. 



In this and the following equations, X denotes the av- 
erage of a random variable X. We eliminated the con- 
tribution of fluctuations of the total atom number n by 
subtracting the conditional expected value of N for a 
given total atom number E(iV|n) — pn, where p denotes 
the emission probability. The latter can be estimated 



from the mean emitted fraction of atoms p w N/n. The 
factor of 2 arises due to the atoms always being emitted 
pairwise. Note, that assuming a Poissonian distribution 
for the shot noise of uncorrelated spontaneous emission 
is generally only valid within the approximation of neg- 
ligible source depletion (N <^ tiq, tiq denoting the pop- 
ulation of the source state), i.e., at early times during 
the emission. Also, we have to neglect fluctuations of 
the relative source population no/n. Those would re- 
sult from fluctuations in temperature and excitation effi- 
ciency, which are inaccessible to independent characteri- 
zation. Hence, an interpretation of the relative variance 
f ^ as being normalized to the shot noise of randomly 
emitted atoms for an ensemble of identical initial condi- 
tions is rather crude. Still it gives the correct order of 
magnitude, and allows for comparison between different 
times during emission, and different excitation strengths. 
In the next step, the contribution of photon shot noise 
in the ffuorescence imaging [26 has to be accounted for. 
This can be done in a similar manner to the calculation 
shown in [5^ (Supplementary Information), and leads to: 



ag. 



Va.r{S-ps)^2^^mS + 2S 
In this equation, s and S denote the total number of pho- 
tons scattered from atoms, and the number of photons 
outside the source region as deflned above, respectively. 
The average number of detected photos emitted by each 
atom is denoted as m, typically m w 12. This means, 
that S = rriN, s = mn, and thus S/s — N/n w p. How- 
ever, Var(S' — ps) > Var(mA^ — pnin). The second term 
represents the photon shot noise of light scattered by the 
atoms, where the factor of 2 arises due to the EMCCD 
camera used for imaging [55^. Background light is taken 
into account by a^ , denoting the variance of the number 
of background photons originating from stray light and 
detector noise within the twin-beam regions. 
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